function vmr = importvmrfromanalyze(img, imeth, hl)
% importvmrfromanalyze  - import a VMR from an Analzye file
%
% FORMAT:       vmr = importvmrfromanalyze(img [, imeth [, histlimit]])
%
% Input fields:
%
%       img         HDR/IMG/NII filename
%       imeth       interpolation method (optional, see flexinterpn)
%       histlimit   limit values on histogram (1x2 double, [0.001 0.999])
%
% Output fields:
%
%       vmr         created VMR object
%
% Note: this function requires the MEX file flexinterpn.

% Version:  v0.8a
% Build:    9102720
% Date:     Oct-27 2009, 8:47 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~ischar(img) || ...
    isempty(img) || ...
    exist(img(:)', 'file') ~= 2
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad input argument.' ...
    );
end
if nargin < 2 || ...
   ~ischar(imeth) || ...
    isempty(imeth)
    imeth = 'cubic';
elseif ~any(strcmpi(imeth(:)', {'cubic', 'lanczos3', 'lanczos8', 'linear', 'nearest'}))
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid interpolation method.' ...
    );
else
    imeth = lower(imeth(:)');
end
if nargin < 3 || ...
   ~isa(hl, 'double')
    hl = [0.001, 0.999];
elseif numel(hl) ~= 2 || ...
    any(isinf(hl) | isnan(hl) | hl < 0 | hl > 4095)
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid histogram limits.' ...
    );
else
    hl = hl(:)';
end
try
    himg = cell(1, 1);
    img = strrep(strrep(img, '.img', '.hdr'), '.IMG', '.HDR');
    himg{1} = BVQXfile(img);
    if numel(himg{1}) ~= 1 ||...
       ~isBVQXfile(himg{1}, 'hdr')
        error('BAD_IMAGECONT');
    end
catch
    clearbvqxobjects(himg);
    error( ...
        'BVQXtools:BadArgument', ...
        'Error loading image.' ...
    );
end

% try to sample volume
try
    vmr16 = himg{1}.SampleBox(struct( ...
        'BBox', [0, 0, 0; 256, 256, 256], 'XYZRes', 1), 1, imeth);
catch
    clearbvqxobjects(himg);
    rethrow(lasterror);
end
mmm = minmaxmean(vmr16, 4);
if mmm(1) ~= 0
    vmr16 = vmr16 - mmm(1);
    mmm = mmm - mmm(1);
end
vmr16 = floor((4095 / mmm(2)) .* vmr16);
vmr16(isinf(vmr16) | isnan(vmr16)) = 0;
vmr16 = uint16(vmr16);

% create object
vmr = bless(BVQXfile('new:vmr'), 1);

% set object properties
updstate = BVQXfile(0, 'updatestate');
BVQXfile(0, 'updatedisable', 'vmr');
vmr.DimX = 256;
vmr.DimY = 256;
vmr.DimZ = 256;
vmr.VMRData = uint8([]);
vmr.VMRData(256, 256, 256) = 0;
vmr.VMRData16 = vmr16;
vmr.FramingCube = 256;
vmr.PosInfoVerified = 1;
vmr.OffsetX = 0;
vmr.OffsetY = 0;
vmr.OffsetZ = 0;
vmr.Slice1CenterX = -127.5;
vmr.Slice1CenterY = 0;
vmr.Slice1CenterZ = 0;
vmr.SliceNCenterX = 127.5;
vmr.SliceNCenterY = 0;
vmr.SliceNCenterZ = 0;
vmr.RowDirX = 0;
vmr.RowDirY = 1;
vmr.RowDirZ = 0;
vmr.ColDirX = 0;
vmr.ColDirY = 0;
vmr.ColDirZ = -1;
vmr.NRows = 256;
vmr.NCols = 256;
vmr.FoVRows = 256;
vmr.FoVCols = 256;
vmr.SliceThickness = 1;
vmr.GapThickness = 0;
vmr.Convention = 1;
vmr.VoxResX = 1;
vmr.VoxResY = 1;
vmr.VoxResZ = 1;
vmr.VoxResInTalairach = 1;
vmr.VoxResVerified = 1;

% reset updatestate flag
BVQXfile(0, 'updatestate', updstate);

% try to perform good limitting
vmr.LimitVMR(struct('recalc8b', true, 'range', hl));
vmr.MeanOriginalValue = floor(mean(vmr.VMRData16(:)));
vmr.MaxOriginalValue = double(max(vmr.VMRData16(:)));
